#include <iostream>
#include<stdio.h>
#include<time.h>
#include <windows.h>
#include <nmmintrin.h>
#include <immintrin.h>//ͷļ
using namespace std;
const int N=50;
void LU(float **m,int n){
    for(int k=0;k<n;k++){
        for(int j=k+1;j<n;j++){
            m[k][j]/=m[k][k];
        }
        m[k][k]=1.0;
        for(int i=k+1;i<n;i++){
                for(int j=k+1;j<n;j++){
                        m[i][j]=m[i][j]-m[i][j]*m[k][j];
        }
        m[i][k]=0;
        }
}
}
void SIMD_SSE(float **m,int N){
 for(int k=0;k<N;k++){
        __m128 vt = _mm_loadu_ps(&m[k][k]);
        for(int j=k+1;j<N-3;j+=4){
            __m128 va= _mm_loadu_ps(&m[k][j]);
            va=_mm_div_ps(va, vt);
            _mm_store_ps(&m[k][j],va);
        }
        for(int j=N-3;j<N;j++)
            m[k][j]=m[k][j]/m[k][k];
        m[k][k]=1.0;
        for(int i=k+1;i<N;i++){
            __m128 vaik=_mm_loadu_ps(&m[i][k]);
             for(int j=k+1;j<N-3;j+=4){
                 __m128 vakj=_mm_loadu_ps(&m[k][j]);
                 __m128 vaij=_mm_loadu_ps(&m[i][j]);
                 __m128 vx=_mm_mul_ps(vakj,vaij);
                 vaij=_mm_sub_ps(vaij,vx);
                 _mm_store_ps(&m[i][j],vaij);
             }
             for(int j=k+1;j<N-3;j+=4){
                 m[i][j]=m[i][j]-m[k][j]*m[i][k];
             }
             m[i][k]=0;
        }
    }
}
int main(){
    float **m=new float*[N];
    for(int i=0;i<N;i++)
      m[i]=new float[N];
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++)
            m[i][j]=rand()%100;
    }

    long long head,tail,freq;
    QueryPerformanceFrequency((LARGE_INTEGER *)&freq);
    QueryPerformanceCounter((LARGE_INTEGER *)&head);
    LU(m,N);
    QueryPerformanceCounter((LARGE_INTEGER *)&tail);
    cout << "BEFORE:"<< (tail-head)* 1000.0 / freq<<"ms"<<endl;

    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++)
            m[i][j]=rand()%100;
    }
    QueryPerformanceCounter((LARGE_INTEGER *)&head);
    SIMD_SSE(m,N);
    QueryPerformanceCounter((LARGE_INTEGER *)&tail );
    cout << "AFTER:"<< (tail-head)* 1000.0 / freq<<"ms"<<endl;

}
